import numpy as np
import matplotlib.pyplot as plt

N = 100000 # Nombre de simulations pour la methode de Monte-Carlo

# Definition des grandeurs mesurees
mCe = 1.012 # Valeur de la masse de sulfate de cerium (en g)
EMT_mCe = 0.0001 # Erreur maximale toleree sur la masse (en g)
VCe = 25.0e-3 # Volume de la fiole jaugee dans laquelle le cerium est dissous (en L)
EMT_VCe = 0.04e-3 # Erreur maximale toleree sur le volume (en L)
MCe = 404.29 # Masse molaire du sulfate de cerium (en g/mol)
Veq = 5.20 # Valeur du volume a l'equivalence (en mL)
EMT_burette = 0.02 # Erreur maximale toleree sur la burette (en mL)
graduation = 0.05 # Demi-graduation de la burette (en mL)
demi_largeur_determination = 0.44/2 # Determination de l'equivalence a l'aide de la demi-largeur a mi-hauteur (en mL)
V0 = 50 # Valeur du volume de solution de fer titree (en mL)
EMT_V0 = 0.05 # Erreur maximale toleree sur le volume de solution de fer titree (en mL)

# Methode de Monte-Carlo
mCe_rd = mCe + np.random.uniform(-EMT_mCe, EMT_mCe, N) + np.random.uniform(-EMT_mCe, EMT_mCe, N)
Veq_rd = Veq + np.random.triangular(-EMT_burette, 0, EMT_burette, N) + np.random.uniform(-graduation, graduation, N) + np.random.uniform(-graduation, graduation, N) + np.random.triangular(-demi_largeur_determination, 0, demi_largeur_determination, N)
VCe_rd = np.random.triangular(VCe-EMT_VCe,VCe, VCe+EMT_VCe, N)
V0_rd = np.random.triangular(V0-EMT_V0, V0, V0+EMT_V0, N)
C = mCe_rd/(MCe*VCe_rd)* Veq_rd/V0_rd
C_calc = mCe/(MCe*VCe)* Veq/V0
s_C = np.std(C)

# Affichage du resultat
print('CM =',format(C_calc,'.3E'),'\u00B1',format(2*s_C,'.1E'),'mol/L a 95 % de confiance')

# Histogramme des differentes valeurs de la concentration en ions ferreux
plt.figure("Tirage")
plt.hist(C, bins = 200, color = 'grey', alpha = 0.7)
plt.xlabel("$C_M$ (mol/L)", fontsize = 20)
plt.ylabel("Occurrences", fontsize = 20)
plt.xticks(fontsize = 20)
plt.yticks(fontsize = 20)
plt.show()

# Diagramme en barres de la contribution relative des incertitudes
C_mCe = mCe_rd/(MCe*VCe)* Veq/V0 # Creation d'un tableau de C avec uniquement mCe qui varie
C_VCe = mCe/(MCe*VCe_rd)* Veq/V0 # Creation d'un tableau de C avec uniquement VCe qui varie
C_Veq = mCe/(MCe*VCe)* Veq_rd/V0 # Creation d'un tableau de C avec uniquement Veq qui varie
C_V0 = mCe/(MCe*VCe)* Veq/V0_rd # Creation d'un tableau de C avec uniquement V0 qui varie
simul= [C_mCe, C_VCe, C_Veq, C_V0]

prop_tot = np.var(C_mCe) + np.var(C_VCe) + np.var(C_Veq) + np.var(C_V0) # Variance totale
valeur =[]
for i in range(len(simul)):
     valeur.append(np.var(simul[i])/prop_tot*100) # Calcul de la contribution relative a la variance

plt.figure("Proportion")
plt.bar([1,2,3,4], valeur, tick_label = ('$m_{Ce}$','$V_{fiole}$' , '$V_{eq}$', '$V_0$'), color = 'grey', alpha = 0.7)
plt.xlabel("Grandeur mesuree", fontsize = 20)
plt.ylabel("Contribution relative de chaque incertitude (%)", fontsize = 20)
for i in range(len(valeur)):
        plt.text(i+1, valeur[i]+1, round(valeur[i],2), ha = 'center', fontsize = 20)
plt.xticks(fontsize = 20)
plt.yticks(fontsize = 20)
plt.show()